Method for propagating pseudo acoustic quasi-p waves in anisotropic media

ABSTRACT

A computer-implemented method for pseudo acoustic quasi-P wave propagation which remain stable in anisotropic media with variable tilt and is not limited to weak anisotropic conditions. The method includes acquiring a seismic exploration volume for a subsurface region of interest, and determining a modeling geometry for the seismic exploration volume. The method further includes propagating at least one wavefield through the seismic exploration volume utilizing the modeling geometry and initial conditions and preventing the accumulation of energy along the axis of symmetry of the seismic exploration volume and ensuring positive stiffness coefficients in the stress-strain relations through the use of finite quasi-S wave velocities thereby producing a stable wavefield. The method includes utilizing the stable wavefield to generate subsurface images of the subsurface region of interest.

FIELD OF THE INVENTION

The present invention relates generally to geophysical prospecting usingseismic signals, and in particular a method for propagatingpseudo-acoustic quasi-P wave propagation in variable tilted anisotropicmedia and using the propagated wavefields for subsurface propertycharacterization.

BACKGROUND OF THE INVENTION

Anisotropy is ubiquitously observed in many oil and gas explorationareas (e.g., the Gulf of Mexico, the North Sea, and offshore WestAfrica) because of preferred ordering of minerals and defects related tostresses. In these regions, often the rock properties can becharacterized as transversely isotropic (“TI”) media with either avertical or tilted axis of symmetry. Wave propagation in anisotropicmedia exhibits different kinematics and dynamics from that in isotropicmedia, thus, it requires anisotropic modeling and migration methods toimage reservoirs properly for oil and gas exploration.

Three-dimensional (“3D”) anisotropic seismic modeling and migration,however, are computationally intensive tasks. Compared to prior artsolutions of full elasticity, modeling and migration based on dispersionrelations are computationally efficient alternatives. In one prior artmethod, Alkhalifah (2000), a pseudo-acoustic approximation for verticaltransversely isotropic (“VTI”) media was introduced. In theapproximation of that prior art method, the phase velocity of shearwaves is set to zero along the vertical axis of symmetry. Thissimplification doesn't eliminate shear waves in other directions asdescribed by Grechka et al. (2004). Based on Alkhalifah's approximation,several space- and time-domain pseudo-acoustic partial differentialequations (PDEs) have been proposed (Alkhalifah, 2000; Zhou et al.,2006; and Du et al., 2008) for seismic modeling and migration in VTImedia. These systems of PDEs are close approximations in kinematics tothe solutions of full elasticity involving vector fields.

As an extension from VTI media, the axis of symmetry of a TI medium canbe tilted (“TTI”) as observed in regions associated with anticlinalstructures and/or thrust sheets. Zhou et al. (2006) extended their VTIpseudo-acoustic equations to a system for 2D TTI media by applying arotation about the axis of symmetry. Consequently the phase velocity ofquasi-SV waves is zero in the direction parallel or perpendicular to thetilted axis. Lesage et al. (2008) further extended Zhou's TTI systemfrom 2D to 3D based on the same phase velocity approximation. However,these prior art pseudo-acoustic modeling and migration methods canbecome numerically unstable due to rapid lateral variations in tiltand/or certain rock properties (when the vertical velocity is greaterthan the horizontal velocity) and result in unstable wave propagation.

As one skilled in the art will appreciate, the plane-wave polarizationvector in isotropic media is either parallel (for P-waves) or orthogonal(for S-waves) to the slowness vector. Except for specific propagationdirections, there are no pure longitudinal and shear waves inanisotropic media. For that reason, in anisotropic wave theory the fastmode is often referred to as the “quasi-P” wave and the slow modes“quasi-S₁” and “quasi-S₂”.

SUMMARY OF THE INVENTION

The present invention provides both a pseudo-acoustic modeling methodand a pseudo-acoustic migration method for anisotropic media. Aspects ofembodiments of the present invention include a computer-implementedmethod for pseudo acoustic quasi-P wave propagation which remain stablein variable-tilt anisotropic media and is not limited to weakanisotropic conditions. The method also includes acquiring a seismicexploration volume for a subsurface region of interest, and determininga modeling geometry for the seismic exploration volume. The methodfurther includes propagating at least one wavefield through the seismicexploration volume utilizing the modeling geometry for initialconditions and preventing the accumulation of energy along the axis ofsymmetry as well as ensuring positive stiffness coefficients in thestress-strain relations through the use of a small finite quasi-S wavevelocity thereby producing a stable wavefield. The method includesutilizing the stable wavefield to generate subsurface images of thesubsurface region of interest.

Another embodiment of the present invention includes a geophysicalseismic migration method comprising the steps of establishing a seismicdata set and a velocity/anisotropy model corresponding to a seismicexploration volume, and for each common shot/receiver record, settingboundary conditions to include excitation from source location(s). Theembodiment further includes propagating wavefields forward according toa pseudo-acoustic wave equation or its equivalents:

$\begin{matrix}{\mspace{79mu} \{ {\begin{matrix}{{\frac{\partial^{2}}{\partial t^{2}}P} = \begin{matrix}{{{v_{P\; 0}^{2}\lbrack {( {{( {1 + {2\; ɛ}} )f_{2}} + f_{1}} ) - {( {f - 1} )f_{3}}} \rbrack}P} +} \\{{v_{P\; 0}^{4}\lbrack {{2\; {f( {\delta - ɛ} )}f_{1}f_{2}} + {( {f - 1} )( {{( {1 + {2\; ɛ}} )f_{2}} + f_{1}} )f_{3}}} \rbrack}Q}\end{matrix}} \\{{\frac{\partial^{2}}{\partial t^{2}}Q} = P}\end{matrix}\mspace{79mu} {where}\text{:}\{ \begin{matrix}{f = {1 - ( \frac{v_{s_{0}}}{v_{p_{0}}} )^{2}}} \\{f_{1} = {{\sin^{2}{\theta_{0}\begin{pmatrix}{{\cos^{2}\varphi_{0}\frac{\partial^{2}}{\partial x^{2}}} +} \\{{\sin^{2}\varphi_{0}\frac{\partial^{2}}{\partial y^{2}}} +} \\{\sin \; 2\varphi_{0}\frac{\partial}{\partial x}\frac{\partial}{\partial y}}\end{pmatrix}}} + {\cos^{2}\theta_{0}\frac{\partial^{2}}{\partial z^{2}}} + {\sin \; 2{\theta_{0}\begin{pmatrix}{{\cos \; \varphi_{0}\frac{\partial}{\partial x}\frac{\partial}{\partial z}} +} \\{\sin \; \varphi_{0}\frac{\partial}{\partial y}\frac{\partial}{\partial z}}\end{pmatrix}}}}} \\{f_{2} = \begin{matrix}{{( {1 - {\sin^{2}\theta_{0}\cos^{2}\varphi_{0}}} )\frac{\partial^{2}}{\partial x^{2}}} + {( {1 - {\sin^{2}\theta_{0}\sin^{2}\varphi_{0}}} )\frac{\partial^{2}}{\partial y^{2}}} +} \\{{\sin^{2}\theta_{0}\frac{\partial^{2}}{\partial z^{2}}} - {\sin^{2}\theta_{0}\sin \; 2\varphi_{0}\frac{\partial}{\partial x}\frac{\partial}{\partial y}} -} \\{\sin \; 2{\theta_{0}( {{\cos \; \varphi_{0}\frac{\partial}{\partial x}\frac{\partial}{\partial z}} + {\sin \; \varphi_{0}\frac{\partial}{\partial y}\frac{\partial}{\partial z}}} )}}\end{matrix}} \\{f_{3} = {{f_{1} + f_{2}} = {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} + \frac{\partial^{2}}{\partial z^{2}}}}}\end{matrix} } } & \lbrack 1\rbrack\end{matrix}$

Vs₀ is the vertical velocity of quasi-SV waves, Vp₀ is the verticalvelocity of quasi-P waves, θ₀ is the tilt of the axis of symmetry withrespect to the vertical in a TI medium, φ₀ is the azimuth of the axis ofsymmetry, ε, δ are the Thomsen anisotropy parameters, P is a scalarwavefield, and Q is an auxiliary function. The embodiment also includesfor each common shot/receiver record, setting boundary conditions toback propagate a recorded shot record, and propagating seismic databackward according to the above pseudo-acoustic wave equations. Theembodiment includes applying imaging conditions such as (but not limitedto) cross correlation between the computed forward wavefields andbackward wavefields or their equivalent Green's functions to derivesubsurface images.

An additional embodiment of the present invention also includes the stepof propagating wavefields or calculating Green's functions by reversetime migration (RTM), Gaussian beam migration, Kirchhoff migration orother wave equation based migrations.

An additional embodiment of the present invention also includes the stepof applying imaging condition involving illumination normalizationand/or reflection-angle domain gather generation and/or phase-amplitudecompensation in addition to cross correlation as options.

An additional embodiment of the present invention also includes the stepof processing common-shot/receiver signals and propagating wavefields inother dependent domains, including but not limited to common offset,common azimuth, and common reflection-angle, and in other modeling andmigration forms, including but not limited to delayed shot, plane-wave,and phase encoding.

An additional embodiment of the present invention also includes the stepof propagating wavefields or calculating Green's functions using otherequivalent terms such as normal moveout velocity, horizontal velocityinstead of Thomsen parameters.

An additional embodiment of the present invention includes a geophysicalseismic migration method comprising the steps of establishing a seismicdata set and a velocity/anisotropy model corresponding to a seismicexploration volume, and for each common shot/receiver record, settingboundary conditions to include excitation from source location(s). Theembodiment also includes propagating wavefields forward according to thefollowing pseudo-acoustic wave equation or its equivalents:

$\begin{matrix}\{ \begin{matrix}{{\frac{\partial^{2}}{\partial t^{2}}P} = \begin{matrix}{{{v_{P\; 0}^{2}\lbrack {( {{( {1 + {2\; ɛ}} )f_{2}} + f_{1}} ) - {( {f - 1} )f_{3}}} \rbrack}P} +} \\{{v_{P\; 0}^{2}\lbrack {{2\; {f( {\delta - ɛ} )}f_{1}f_{2}} + {( {f - 1} )( {{( {1 + {2\; ɛ}} )f_{2}} + f_{1}} )f_{3}}} \rbrack}Q}\end{matrix}} \\{{\frac{\partial^{2}}{\partial t^{2}}Q} = {v_{p_{0}}^{2}P}}\end{matrix}  & \lbrack 2\rbrack\end{matrix}$

The embodiment further includes for each common shot/receiver record,setting boundary conditions to back propagate a recorded shot record,and propagating seismic data backward according to the abovepseudo-acoustic wave equations. The embodiment includes applying imagingconditions such as (but not limited to) cross correlation between thecomputed forward wavefields and backward wavefields or their equivalentGreen's functions to derive subsurface images.

Different embodiments of the present invention may utilize otherpseudo-acoustic wave equations to propagate wavefields forward ingeophysical seismic migration.

For example, one embodiment of the present invention includespropagating wavefields forward according to the pseudo-acoustic waveequation below or its equivalents:

$\begin{matrix}\{ \begin{matrix}{{\omega^{2}P} = \begin{matrix}{{{v_{P\; 0}^{2}\lbrack {( {{( {1 + {2\; ɛ}} )f_{2}} + f_{1}} ) - {( {f - 1} )f_{3}}} \rbrack}P} +} \\{{v_{P\; 0}^{4}\lbrack {{2\; {f( {\delta - ɛ} )}f_{1}f_{2}} + {( {f - 1} )( {{( {1 + {2\; ɛ}} )f_{2}} + f_{1}} )f_{3}}} \rbrack}Q}\end{matrix}} \\{{\omega^{2}Q} = P}\end{matrix}  & \lbrack 3\rbrack\end{matrix}$

where ω is the angular frequency.

A further embodiment of the present invention that is utilized forgeophysical seismic migration includes propagating wavefields forwardaccording to the pseudo-acoustic wave equation or its equivalents:

$\begin{matrix}\{ \begin{matrix}{{\frac{\partial^{2}}{\partial t^{2}}P} = \begin{matrix}{{{v_{P\; 0}^{2}\lbrack {( {{( {1 + {2\; ɛ}} )f_{2}} + f_{1}} ) - {( {f - 1} )f_{3}}} \rbrack}P} +} \\{{v_{P\; 0}^{4}{f\begin{bmatrix}{( {1 + {2\; \delta}} ) -} \\( {1 + {2ɛ}} )\end{bmatrix}}f_{1}Q} + {{{v_{P\; 0}^{4}( {f - 1} )}\begin{bmatrix}( {1 + {2\; ɛ}} ) \\{f_{2} + f_{1}}\end{bmatrix}}R}}\end{matrix}} \\{{\frac{\partial^{2}}{\partial t^{2}}Q} = {f_{2}P}} \\{{\frac{\partial^{2}}{\partial t^{2}}R} = {f_{3}P}}\end{matrix}  & \lbrack 4\rbrack\end{matrix}$

where Q and R are auxiliary functions.

Another embodiment of the present invention includes a geophysicalseismic migration method comprising the steps of establishing a seismicdata set and a velocity/anisotropy model corresponding to a seismicexploration volume, and for each common shot/receiver record, settingboundary conditions to include excitation from source location(s). Theembodiment also includes propagating wavefields forward according to apseudo-acoustic wave equation and its equivalent formulations for tiltedmedia:

$\begin{matrix}\{ \begin{matrix}{{\frac{\partial}{\partial t}P} = {{( {1 + {2\; \eta}} ){v_{nmo}^{2}( {{\frac{\partial}{\partial x}U} + {\frac{\partial}{\partial y}V}} )}} + {\frac{v_{p_{0}}^{2}}{1 + {2\eta}}\frac{\partial}{\partial z}R}}} \\{{\frac{\partial}{\partial t}U} = {\frac{\partial}{\partial x}P}} \\{{\frac{\partial}{\partial t}V} = {\frac{\partial}{\partial y}P}} \\{{\frac{\partial}{\partial t}Q} = {\frac{v_{p_{0}}^{2}}{1 + {2\eta}}\frac{\partial}{\partial z}R}} \\{{\frac{\partial}{\partial t}R} = {{\frac{\partial}{\partial z}P} + {2\eta \frac{\partial}{\partial z}Q}}}\end{matrix}  & \lbrack 5\rbrack\end{matrix}$

v_(p) ₀ is the vertical velocity of quasi-P waves, v_(nmo)=v_(p) ₀√{square root over (1+2δ)} is the normal-moveout velocity of quasi-Pwaves, η=(ε−δ)/(1+2δ) is the Alkhalifah-Tsvankin anisotropy parameter(expressed in terms of the Thomsen anisotropy parameters ε and δ), P isa scalar wavefield, and U, V, Q, and R are auxiliary functions. Theembodiment further includes for each common shot/receiver record,setting boundary conditions to back propagate a recorded shot record,and propagating seismic data backward according to the abovepseudo-acoustic wave equations. The embodiment includes applying imagingconditions such as cross correlation between the computed forward andbackward wavefields or their equivalent Green's functions to derivesubsurface images.

Different embodiments of the present invention for geophysical seismicmigration may utilize other pseudo-acoustic wave equations to propagatewavefields forward for tilted media. For example, one embodiment of thepresent invention includes propagating wavefields forward according to apseudo-acoustic wave equation and its equivalent formulations for tiltedmedia:

$\begin{matrix}\{ {\begin{matrix}{{\frac{\partial^{2}}{\partial t^{2}}P} = \begin{matrix}{{\begin{bmatrix}{{( {1 + {2\eta}} )v_{nmo}^{2}} +} \\{av}_{po}^{2}\end{bmatrix}g_{2}P} + {( {1 + a} )v_{po}^{2}g_{1}P} -} \\{{{a( {1 + {2\eta}} )}v_{nmo}^{2}g_{2}Q} - {\begin{bmatrix}{{( {{2\eta} + a} )v_{nmo}^{2}} +} \\{av}_{po}^{2}\end{bmatrix}g_{1}Q} - {{av}_{po}^{2}g_{1}R}}\end{matrix}} \\{{\frac{\partial^{2}}{\partial t^{2}}Q} = {v_{po}^{2}g_{2}P}} \\{{\frac{\partial^{2}}{\partial t^{2}}R} = {v_{po}^{2}g_{1}P}}\end{matrix}\mspace{79mu} {where}\mspace{79mu} \{ \begin{matrix}{a = {{1 - f} = ( {v_{s_{0}}/v_{p_{0}}} )^{2}}} \\{g_{1} = \frac{\partial^{2}}{\partial z^{2}}} \\{g_{1} = {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}}}\end{matrix} }  & \lbrack 6\rbrack\end{matrix}$

v_(p) ₀ is the vertical velocity of quasi-P waves, v_(nmo)=v_(p) ₀√{square root over (1+2δ)} is the normal-moveout velocity of quasi-Pwaves, a is the square of the shear-wave to P-wave velocity ratio,η=(ε−δ)/(1+2δ) is the Alkhalifah-Tsvankin anisotropy parameter(expressed in terms of the Thomsen anisotropy parameters ε and δ), and Qand R are auxiliary functions.

A further embodiment of the present invention that is utilized forgeophysical seismic migration includes propagating wavefields forwardaccording to a pseudo-acoustic wave equation and its equivalentformulation for tilted media or their derivativeformulations/equivalents:

$\begin{matrix}{{{\frac{\partial^{4}}{\partial t^{4}}F} - \{ {{\lbrack {{( {1 + {2\eta}} )v_{nmo}^{\;^{2}}} + {av}_{po}^{2}} \rbrack ( {{\frac{\partial^{4}}{{\partial x^{2}}{\partial t^{2}}}F} + {\frac{\partial^{4}}{{\partial y^{2}}{\partial t^{2}}}F}} )} + {( {1 + a} ){v_{po}^{2}( {\frac{\partial^{4}}{{\partial z^{2}}{\partial t^{2}}}F} )}}} \} + {{a( {1 + {2\eta}} )}v_{nmo}^{2}{v_{po}^{2}( {{\frac{\partial^{4}}{\partial x^{4}}F} + {\frac{\partial^{4}}{{\partial x^{2}}{\partial y^{2}}}F} + {\frac{\partial^{4}}{\partial y^{4}}F}} )}} + {\lbrack {{( {{2\eta} + a} )v_{nmo}^{2}v_{po}^{2}} + {av}_{po}^{2}} \rbrack ( {{\frac{\partial^{4}}{{\partial x^{2}}{\partial z^{2}}}F} + {\frac{\partial^{4}}{{\partial y^{2}}{\partial z^{2}}}F}} )} + {{av}_{po}^{4}( {\frac{\partial^{4}}{\partial z^{2}}F} )}} = 0} & \lbrack 7\rbrack\end{matrix}$

where F is a scalar wavefield.

Another embodiment of the present invention includes a geophysicalseismic migration method comprising the steps of establishing a seismicdata set and a velocity/anisotropy model corresponding to a seismicexploration volume, and for each common shot/receiver record, settingboundary conditions to include excitation from source location(s). Theembodiment also includes propagating wavefields forward according to apseudo-acoustic wave equation or its derivativeformulations/equivalents:

$\begin{matrix}{{\frac{\partial^{4}P}{\partial t^{4}} - {{v_{po}^{2}\lbrack {{( {1 + {2ɛ}} )f_{2}} + f_{1} - {( {f - 1} )f_{3}}} \rbrack}\frac{\partial^{2}P}{\partial t^{2}}} + {{v_{po}^{4}\lbrack {{2\; {f( {ɛ - \delta} )}f_{1}f_{2}} - {( {f - 1} )( {{( {1 + {2ɛ}} )f_{2}} + f_{1}} )f_{3}}} \rbrack}P}} = 0} & \lbrack 8\rbrack\end{matrix}$

The embodiment further includes for each common shot/receiver record,setting boundary conditions to back propagate a recorded shot record,and propagating seismic data backward according to the abovepseudo-acoustic wave equations. The embodiment includes applying imagingconditions such as (but not limited to) cross correlation between thecomputed forward and backward wavefields or their equivalent Green'sfunctions to derive subsurface images.

Another embodiment of the present invention includes a geophysicalseismic modeling method comprising the steps of establishing avelocity/anisotropy model corresponding to a seismic exploration volume,and for each shot, setting initial conditions of wavefields. Theembodiment also includes propagating wavefields forward according to apseudo-acoustic wave equation or its equivalents:

$\begin{matrix}\{ \begin{matrix}{{\frac{\partial^{2}}{\partial t^{2}}P} = \begin{matrix}{{{v_{P\; 0}^{2}\lbrack {( {{( {1 + {2\; ɛ}} )f_{2}} + f_{1}} ) - {( {f - 1} )f_{3}}} \rbrack}P} +} \\{{{v_{P\; 0}^{4}\lbrack {{2\; {f( {\delta - ɛ} )}f_{1}f_{2}} + {( {f - 1} )( {{( {1 + {2\; ɛ}} )f_{2}} + f_{1}} )f_{3}}} \rbrack}Q} +} \\{{\delta ( {\overset{arrow}{x} - {\overset{arrow}{x}}_{s}} )}{w(t)}}\end{matrix}} \\{{\frac{\partial^{2}}{\partial t^{2}}Q} = {P + {{\delta ( {\overset{arrow}{x} - {\overset{arrow}{x}}_{s}} )}{w(t)}}}}\end{matrix}  & \lbrack 9\rbrack\end{matrix}$

where w(t) is a source function, and {right arrow over (x)}_(s) is thevector of the source location. The source term and its form of insertioncan be changed without affecting the governing PDEs.

Another embodiment of the present invention that is utilized forgeophysical seismic modeling includes propagating wavefields forwardaccording to a pseudo-acoustic wave equation (equation 5) and itsequivalent formulations for tilted media.

Another embodiment of the present invention that is utilized forgeophysical seismic modeling includes propagating wavefields forwardaccording to a pseudo-acoustic wave equation (equation 6) and itsequivalent formulations for tilted media.

It should also be appreciated that the present invention is intended tobe used with a system which includes, in general, an electronicconfiguration including at least one processor, at least one memorydevice for storing program code or other data, an optional video monitoror other display device (i.e., a liquid crystal display) and at leastone input device. The processor is preferably a microprocessor ormicrocontroller-based platform which is capable of displaying images andprocessing complex mathematical algorithms The memory device can includerandom access memory (RAM) for storing event or other data generated orused during a particular process associated with the present invention.The memory device can also include read only memory (ROM) for storingthe program code for the controls and processes of the presentinvention.

One such embodiment includes a system configured to perform pseudoacoustic quasi-P wave propagation which remain stable in variable tiltanisotropic media and is not limited to weak anisotropic conditions. Thesystem includes a data storage device having computer readable dataincluding a seismic exploration volume for a subsurface region ofinterest, and a processor, configured and arranged to execute machineexecutable instructions stored in a processor accessible memory forperforming a method. The method for this particular embodiment includesdetermining a modeling geometry for the seismic exploration volume, andpropagating at least one wavefield through the seismic explorationvolume utilizing the modeling geometry for initial conditions andpreventing the accumulation of energy along the axis of symmetry ofanisotropic regions within the seismic exploration volume and ensuringpositive stiffness coefficients in the stress-strain relations therebyproducing a stable wavefield. The method further includes utilizing thestable wavefield to generate subsurface images of the subsurface regionof interest.

These and other objects, features, and characteristics of the presentinvention, as well as the methods of operation and functions of therelated elements of structure and the combination of parts and economiesof manufacture, will become more apparent upon consideration of thefollowing description and the appended claims with reference to theaccompanying drawings, all of which form a part of this specification,wherein like reference numerals designate corresponding parts in thevarious Figures. It is to be expressly understood, however, that thedrawings are for the purpose of illustration and description only andare not intended as a definition of the limits of the invention. As usedin the specification and in the claims, the singular form of “a”, “an”,and “the” include plural referents unless the context clearly dictatesotherwise.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow chart illustrating a method in accordance with one ormore embodiments of the present invention.

FIG. 2 is a flow chart illustrating a method in accordance with one ormore embodiments of the present invention.

FIG. 3 is a flow chart illustrating a method in accordance with one ormore embodiments of the present invention.

FIG. 4 illustrates exemplary wave propagation modeling according to theprior art, Alkhalifah's approximation where f=1.

FIG. 5 illustrates exemplary wave propagation modeling according to oneembodiment of the present invention.

FIG. 6 illustrates exemplary wave propagation modeling according to oneembodiment of the present invention where Vs₀/Vp₀=0.01.

FIG. 7 illustrates exemplary wave propagation modeling according to theprior art, Alkhalifah's approximation where Vs₀=0.

FIG. 8 illustrates an exemplary phase velocity distribution according tothe prior art, Alkhalifah's approximation.

FIG. 9 illustrates an exemplary group velocity distribution according tothe prior art, Alkhalifah's approximation.

FIG. 10 illustrates an exemplary phase velocity distribution for oneembodiment of the present invention where Vs₀/Vp₀=0.01.

FIG. 11 illustrates an exemplary group velocity distribution for oneembodiment of the present invention where Vs₀/Vp₀=0.01.

FIG. 12 illustrates an exemplary wave propagation modeling in a mediumwith a variable tilted axis of symmetry, according to one embodiment ofthe present invention utilizing a first-order 5×5 PDE system.

FIG. 13 illustrates a schematic diagram of the geometry that is used inone embodiment of the present invention.

FIG. 14 illustrates is a schematic illustration of an embodiment of asystem for performing methods in accordance with embodiments of thepresent invention.

DETAILED DESCRIPTION

One embodiment of the present invention is illustrated in FIG. 1,wherein a flow chart 10 describes a method for propagating quasi-P waveswhich remain stable in anisotropic media with variable tilt. The presentinvention is not limited to weak anisotropic conditions. This particularembodiment includes acquiring a seismic exploration volume of asubsurface region of interest 12, and determining a modeling geometryfor the seismic exploration volume 14. The embodiment further includespropagating at least one wavefield through the seismic explorationvolume utilizing the modeling geometry for initial conditions andpreventing the accumulation of energy along the axis of symmetry for theseismic exploration volume and ensuring positive stiffness coefficientsin the stress-and-strain relations utilizing finite quasi-S wavevelocities thereby producing a stable wavefield 16. The stable wavefieldcan then be utilized to generate subsurface images of the subsurfaceregion of interest 18.

As one in skilled in the art will appreciate, differing embodiments ofthe present invention may provide a pseudo-acoustic modeling method or apseudo-acoustic migration method for anisotropic media. For example,FIG. 2 illustrates a flowchart 20 for one embodiment of apseudo-acoustic modeling method for wave propagation in anisotropicmedia with variable tilt, wherein the method is not limited to weakanisotropic conditions. That embodiment includes acquiring a seismicexploration volume for a subsurface region of interest 22 anddetermining a modeling geometry for the seismic exploration volume 24.The embodiment also includes propagating at least one wavefield throughthe seismic exploration volume utilizing the modeling geometry forinitial conditions, wherein the artificial quasi-shear wave velocity isgreater or equal to zero along the axis of symmetry for the seismicexploration volume thereby preventing the accumulation of energy alongthe axis of symmetry thereby producing a stable wavefield 26. The stablewavefield can then be utilized to generate subsurface images for thesubsurface region of interest 28.

FIG. 3 illustrates a flowchart 30 for another embodiment of the presentinvention that can be used for pseudo-acoustic migration. Thatembodiment includes acquiring a seismic exploration volume for asubsurface region of interest 32 and determining a model geometry forthe seismic exploration volume 34. The embodiment also includespropagating at least one wavefield through the seismic explorationvolume utilizing the modeling geometry for initial conditions, whereinquasi-shear wave energy does not accumulate along the axis of symmetryfor the seismic exploration volume thereby producing a stable wavefield36. The stable wavefield can then be utilized to generate subsurfaceimages for the subsurface region of interest 38.

The present invention provides several advantages relative toconventional acoustic, anisotropic modeling and migration. The presentinvention provides a stable way of wave propagation in TI media withvariable tilt, thus simulated wavefield propagation and images ofreflectivity can be obtained. Prior art pseudo-acoustic modeling andmigration methods are based on Alkhalifah's approximation in which thephase velocity of shear waves is set to zero along the axis of symmetry.Although the prior art methods can work in a constant-tilt TI medium,the zero-speed shear-waves can make wave propagation unstable (i.e.amplitudes become unbounded) in areas where tilt variations can locallyconcentrate high energy near the axis of symmetry. FIG. 4 shows thatprior art methods (f=1) are unstable 40 in a variable-tilt medium (e.g.near the crest of an anticlinal structure). On the contrary, FIG. 5shows that wave propagation based on the present invention (f=0.98)remains stable 42 in the same medium. In addition, the present inventioncan provide the flexibility of controlling shear- to P-wave velocityratios to optimize the results of modeling and migration. For example,shear- and P-wave velocity ratios can be set close to the actual valuesto approximate the kinematics in elastic wave propagation. Furthermore,in certain rocks, the vertical velocity may be greater than thehorizontal velocity with respect to the axis of symmetry. In such acase, wave equations based on Alkhalifah's approximation will result innegative stiffness matrix thereby producing unstable wavefieldsregardless numerical implementation algorithms. The present inventioncan use a finite shear-wave velocity to ensure positive stiffnesscoefficients in the stress-strain relations thereby generating stablewave propagation.

In prior art methods based on Alkhalifah's approximation, equatingshear-wave phase velocity to zero does not eliminate shear waves.Instead, high energy concentrates near the axis of symmetry. The onlyexception is elliptic anisotropy, (i.e., ε=δ), for which shear wavesvanish everywhere. In the present invention, vertical shear velocity isrelaxed from being zero, hence, the energy is less concentrated near theaxis of symmetry. Shear waves will not vanish because of the presence ofadditional cross derivatives even if the conditions of ellipticalanisotropy are satisfied.

In terms of computational cost, the PDEs utilized by embodiments of thepresent invention involve additional spatial derivative terms to becomputed compared to prior art methods. In areas with variable tilt, theadditional workload associated with non-zero Vs₀ is necessary to achievestability and reliability required by seismic modeling and migration. Inareas of nearly constant or very smooth tilt, the additional workloadmay be skipped.

It will be clear to one skilled in the art that the above embodimentsmay be altered in many ways without departing from the scope of theinvention. For example, as is apparent to one skilled in the art,different initial conditions or boundary conditions or a differentlinear combination of the PDEs in the present invention can be used inmodeling and migration as convenient.

In one embodiment of the present invention, the anisotropic modelingmethod includes establishing a velocity and anisotropy modelcorresponding to a seismic exploration volume; setting initialconditions such as source excitation; propagating waves in transverselyisotropic media with a tiled or vertical axis of symmetry, according toeq. [1] or its equivalent. For forward modeling, a source function ofthe form δ({right arrow over (x)}−{right arrow over (x)}_(s))w(t) needsto be introduced in the right side of equations in eq. [1] or eq. [2],where {right arrow over (x)}_(s) is the source location, and w(t) is asource wavelet.

In the above-described embodiment of the present invention, the verticalshear-wave velocity in eq. [3] can be non-zero (therefore f can bedifferent from 1) in contrast to the prior art method approximationwhere f rounds off to 1. Accordingly, the phase velocity of shear wavesin the direction of both parallel and perpendicular to the axis ofsymmetry can be non-zero in the present invention. In a medium withvariable tilt, the finite speed of quasi-shear waves can avoid localconcentration of high energy which often occurs in the vicinity of theaxis of symmetry. The present invention does not require weak anisotropyassumptions.

Utilizing the above PDEs, other embodiments of the present invention canbe derived for anisotropic media. If tilt θ₀=0, the above PDEs simplifyto a 3D VTI system, and similarly a 3D HTI system when φ₀=0. Asecond-order 3×3 system for 3D VTI media can take of the form inequation 6. This system of PDEs is extendable to its equivalentformulation for a tilted TI medium by replacing g₁ and g₂ given in eq.[6] by f₁ and f₂ given in eq. [1].

As an alternative to the PDEs in eq. [1] or [2] or [4] when f=1, thefirst-order 5×5 system of PDEs in eq. [5] is hyperbolic and stable in aTI medium with variable tilt. This embodiment of the present inventionis symmetrizably hyperbolic (well-posed, even with variablecoefficients). This system is also extendable to variable-tilt TTI. Theabove complete first-order 5×5 system of PDEs in 3D reduces to 4×4 in2D.

As described above, additional embodiments of the present invention alsoprovide pseudo-acoustic migration methods. One embodiment includes thesteps of: establishing a seismic data set and a velocity/anisotropymodel corresponding to a seismic exploration volume; setting boundaryconditions of wave propagation; propagating waves from source excitationand recorded seismic data separately in anisotropic media according toeq. [1], eq. [2], eq. [4], or eq. [6], or their equivalents; andapplying imaging conditions such as, but not limited to, crosscorrelation between the two propagated wavefields to obtain subsurfaceimages. Different initial and/or boundary conditions can be appliedwithout affecting the scope of this invention. An exemplary boundarycondition (e.g., based on eq. [1]) for propagating a source wavelet isas follows:

$\begin{matrix}\{ \begin{matrix}{{P( {x,y,{{z = 0};t}} )} = {{\delta ( {\overset{arrow}{x} - {\overset{arrow}{x}}_{s}} )}{\int_{0}^{t}{{w( t^{\prime} )}\ {t^{\prime}}}}}} \\{{Q( {x,y,{{z = 0};t}} )} = {{\delta ( {\overset{arrow}{x} - {\overset{arrow}{x}}_{s}} )}{\int_{0}^{t}{{w( t^{\prime} )}\ {t^{\prime}}}}}}\end{matrix}  & \lbrack 10\rbrack\end{matrix}$

and the boundary condition for reverse time extrapolation of seismicdata is as follows:

$\begin{matrix}\{ \begin{matrix}{{P( {x,y,{{z = 0};t}} )} = {D( {x,y,x_{s},{y_{s};t}} )}} \\{{Q( {x,y,{{z = 0};t}} )} = {D( {x,y,x_{s},{y_{s};t}} )}}\end{matrix}  & \lbrack 11\rbrack\end{matrix}$

where w(t) is a source function, x_(s) is the location of source,D(x,y,x_(s),y_(s);t) is a shot record to migrate.

The following example illustrates a further embodiment of the presentinvention:

1. Establishing Fourth-Order Dispersion Relations for Quasi-P Wave inVTI Media

Tsvankin's phase velocity relations for VTI media for which v_(s0) isnot set to zero lead to the dispersion relation:

$\begin{matrix}{{\omega^{4} - {B\; \omega^{2}} + C} = 0} & \lbrack 12\rbrack \\{{where}\text{:}} & \; \\\{ \begin{matrix}{B = {{\lbrack {{( {1 + {2\eta}} )v_{nmo}^{2}} + {a\; v_{p_{0}}^{2}}} \rbrack k_{h}^{2}} + {( {1 + a} )v_{p_{0}}^{2}k_{z}^{2}}}} \\{C = {{{a( {1 + {2\eta}} )}v_{nmo}^{2}v_{p_{0}}^{2}k_{h}^{4}} + {\lbrack {{( {{2\eta} + a} )v_{nmo}^{2}v_{p_{0}}^{2}} + {a\; v_{p_{0}}^{4}}} \rbrack k_{h}^{2}k_{z}^{2}} + {a\; v_{p_{0}}^{4}k_{z}^{4}}}}\end{matrix}  & \lbrack 13\rbrack\end{matrix}$

ω is angular frequency, k_(z) is the vertical wavenumber, and k_(h)²=k_(x) ²+k_(y) ² is the [square of the] magnitude of thehorizontalwavenumber vector (k_(x), k_(y)). Eq. [12] admits two pairs ofsolutions:

$\begin{matrix}\{ \begin{matrix}{\omega_{{q\; P} \pm} = {\pm \sqrt{\frac{B + \sqrt{B^{2} - {4C}}}{2}}}} \\{\omega_{{q\; S} \pm} = {\pm \sqrt{\frac{B - \sqrt{B^{2} - {4C}}}{2}}}}\end{matrix}  & \lbrack 14\rbrack\end{matrix}$

ω_(qP±) correspond to quasi-P waves; ω_(qS±) correspond to quasi-SVwaves.

2. Establishing a Fourth-Order PDE for Quasi-P Wave in VTI Media

Applying eq. [12] to the wavefield {circumflex over(F)}(k_(x),k_(x),k_(x),ω) in the Fourier domain and taking the inverseFourier transform (F(x,y,z,t)) provide:

$\begin{matrix}{{\frac{\partial^{4}}{\partial t^{4}}F} - \{ {{\lbrack {{( {1 + {2\eta}} )v_{nmo}^{2}} + {a\; v_{p_{0}}^{2}}} \rbrack ( {{\frac{\partial^{4}}{{\partial x^{2}}{\partial t^{2}}}F} + {\frac{\partial^{4}}{{\partial y^{2}}{\partial t^{2}}}F}} )} + {( {1 + a} ){v_{p_{0}}^{2}( {\frac{\partial^{4}}{{\partial z^{2}}{\partial t^{2}}}F} )}}} \} + {{a( {1 + {2\eta}} )}v_{nmo}^{2}{v_{p_{0}}^{2}( {{\frac{\partial^{4}}{\partial x^{4}}F} + {\frac{\partial^{4}}{{\partial x^{2}}{\partial y^{2}}}F} + {\frac{\partial^{4}}{\partial y^{4}}F}} )}} + {\quad{{{\lbrack {{( {{2\eta} + a} )v_{nmo}^{2}v_{p\; 0}^{2}} + {a\; v_{p_{0}}^{2}}} \rbrack ( {{\frac{\partial^{4}}{{\partial x^{2}}{\partial z^{2}}}F} + {\frac{\partial^{4}}{{\partial y^{2}}{\partial z^{2}}}F}} )} + {a\; {v_{p_{0}}^{4}( {\frac{\partial^{4}}{\partial z^{2}}F} )}}} = 0}}} & \lbrack 15\rbrack\end{matrix}$

3. Establishing a Second-Order 3×3 System of PDEs for VTI Media

Let:

$\begin{matrix}{{P( {x,y,z,t} )} = {\frac{\partial^{2}}{\partial t^{2}}{F( {x,y,z,t} )}}} & \lbrack 16\rbrack\end{matrix}$

where F(x,y,z,t) is a wavefield satisfying eq. [15]. Assuming theinitial conditions

$\begin{matrix}{{F( {x,y,z,{t = 0}} )} \equiv {\frac{\partial}{\partial t}{F( {x,y,z,{t = 0}} )}} \equiv 0} & \lbrack 17\rbrack\end{matrix}$

leads to:

F(x,y,z,t)=∫₀ ^(t)∫₀ ^(t′) P(x,y,z,t″)dt″dt′

Let:

$\begin{matrix}\{ \begin{matrix}\begin{matrix}{{Q( {x,y,z,t} )} = {v_{p_{0}}^{2}( {{\frac{\partial^{2}}{\partial x^{2}}F} + {\frac{\partial^{2}}{\partial y^{2}}F}} )}} \\{= {\int_{0}^{t}{\int_{0}^{t^{\prime}}{( {{\frac{\partial^{2}}{\partial x^{2}}P} + {\frac{\partial^{2}}{\partial y^{2}}P}} )( {x,y,z,t^{''}} ){t^{''}}{t^{\prime}}}}}}\end{matrix} \\\begin{matrix}{{R( {x,y,z,t} )} = {v_{p_{0}}^{2}( {\frac{\partial^{2}}{\partial z^{2}}F} )}} \\{= {\int_{0}^{t}{\int_{0}^{t^{\prime}}{( {\frac{\partial^{2}}{\partial z^{2}}P} )( {x,y,z,t^{''}} ){t^{''}}{t^{\prime}}}}}}\end{matrix}\end{matrix}  & \lbrack 19\rbrack\end{matrix}$

Eq. [15] is then equivalent to the second-order 3×3 system of PDEs byeq. [4].

FIG. 6 shows a wavefront propagation in a VTI medium using the abovePDEs for that particular embodiment of the present invention. Comparedto the wavefronts based on the prior art methods (illustrated in FIG.7), the outer qP-wavefront (44 in FIG. 6 and 48 in FIG. 7) remainsalmost identical, but the inner qSV-wavefront (46 in FIG. 6 and 50 inFIG. 7) has a different form from a diamond shape. FIG. 8 and FIG. 9show the phase and group velocities, respectively, according toAlkhalifah's (prior art) approximation. In contrast, FIG. 10 and FIG. 11show the phase and group velocities, respectively, according to anembodiment of the present invention. Compared to Alkhalifah'sapproximation, the phase velocities of qSV waves are relaxed from beingzero along the axis of symmetry. Consequently, the maximum values ofgroup velocities or high energy are not so focused along the axis ofsymmetry as in the prior art methods. The same observations areapplicable to a constant-tilt TTI medium by applying a rotation aboutthe tilt.

4. Establishing a First-Order 5×5 System of PDEs for VTI Media

Du et al. (2008) presents the following second-order 2×2 system of PDEsfor VTI media with v_(s0)=0:

$\begin{matrix}\{ \begin{matrix}{{\frac{\partial}{\partial t}p} = {{( {1 + {2\eta}} )v_{nmo}^{2}g_{2}p} + {v_{p_{0}}^{2}g_{1}q}}} \\{{\frac{\partial}{\partial t}q} = {{v_{nmo}^{2}g_{2}p} + {v_{p_{0}}^{2}g_{1}q}}}\end{matrix}  & \lbrack 20\rbrack\end{matrix}$

where g₁ and g₂ are given after eq. [6], p is a scalar wavefield, and qis an auxiliary function. A new wavefield P is defined and new auxiliaryfunctions U, V, Q, and R by:

$\begin{matrix}{{P = {\frac{\partial}{\partial t}p}},\mspace{14mu} {U = {\frac{\partial}{\partial x}p}},\mspace{14mu} {V = {\frac{\partial}{\partial y}p}},\mspace{14mu} {Q = {\frac{\partial}{\partial t}( \frac{{( {1 + {2\eta}} )q} - p}{2\eta} )}},\mspace{14mu} {R = {( {1 + {2\eta}} )\frac{\partial}{\partial z}{q.}}}} & \lbrack 21\rbrack\end{matrix}$

Then equation 5 is a complete first-order 5×5 system of PDEs. Thissystem can be shown to be hyperbolic by symmetrizing it. Let:

$\begin{matrix}{{\overset{\sim}{P} = P},\mspace{14mu} {\overset{\sim}{U} = {v_{nmo}\sqrt{1 + {2\eta}}U}},\mspace{14mu} {\overset{\sim}{V} = {v_{nmo}\sqrt{1 + {2\eta}}V}},\mspace{14mu} {\overset{\sim}{Q} = {\sqrt{2\eta}Q}},\mspace{14mu} {\overset{\sim}{R} = {\frac{v_{p_{0}}}{\sqrt{1 + {2\eta}}}{R.}}}} & \lbrack 22\rbrack \\{\mspace{20mu} {{Then}\text{:}}} & \; \\{\mspace{20mu} {{\frac{\partial}{\partial t}\begin{bmatrix}\overset{\sim}{P} \\\overset{\sim}{U} \\\overset{\sim}{V} \\\overset{\sim}{Q} \\\overset{\sim}{R}\end{bmatrix}} = {{M_{x}{\frac{\partial}{\partial x}\begin{bmatrix}\overset{\sim}{P} \\\overset{\sim}{U} \\\overset{\sim}{V} \\\overset{\sim}{Q} \\\overset{\sim}{R}\end{bmatrix}}} + {M_{y}{\frac{\partial}{\partial y}\begin{bmatrix}\overset{\sim}{P} \\\overset{\sim}{U} \\\overset{\sim}{V} \\\overset{\sim}{Q} \\\overset{\sim}{R}\end{bmatrix}}} + {M_{z}{\frac{\partial}{\partial z}\begin{bmatrix}\overset{\sim}{P} \\\overset{\sim}{U} \\\overset{\sim}{V} \\\overset{\sim}{Q} \\\overset{\sim}{R}\end{bmatrix}}}}}} & \lbrack 23\rbrack \\{\mspace{20mu} {{where}\text{:}}} & \; \\{\mspace{20mu} {M_{x} = \begin{bmatrix}0 & {v_{nmo}\sqrt{1 + {2\eta}}} & 0 & 0 & 0 \\{v_{nmo}\sqrt{1 + {2\eta}}} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0\end{bmatrix}}} & \lbrack 24\rbrack \\{\mspace{20mu} {M_{y} = \begin{bmatrix}0 & 0 & {v_{nmo}\sqrt{1 + {2\eta}}} & 0 & 0 \\0 & 0 & 0 & 0 & 0 \\{v_{nmo}\sqrt{1 + {2\eta}}} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0\end{bmatrix}}} & \lbrack 25\rbrack \\{\mspace{20mu} {M_{z} = \begin{bmatrix}0 & 0 & 0 & \frac{v_{p_{0}}}{\sqrt{1 + {2\eta}}} & 0 \\0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 \\\frac{v_{p_{0}}}{\sqrt{1 + {2\eta}}} & 0 & 0 & 0 & \frac{v_{p_{0}}\sqrt{2\eta}}{\sqrt{1 + {2\eta}}} \\0 & 0 & 0 & \frac{v_{p_{0}}\sqrt{2\eta}}{\sqrt{1 + {2\eta}}} & 0\end{bmatrix}}} & \lbrack 26\rbrack\end{matrix}$

In 2D, the variable V is eliminated and the third equation is deleted,resulting in a first-order 4×4 system.

FIG. 12 shows stable wavefront propagation governed by such first-orderPDEs in a variable-tilt medium.

5. Establishing Fourth-Order Dispersion Relations in TTI Media

By relaxing Alkhalifah's approximation that Vs₀=0 (or f=1) along theaxis of symmetry, the following equation can be derived from Tsvankin'sphase velocity relations (2001):

$\quad\{ \begin{matrix}{\frac{\cos^{2}\overset{\_}{\theta}}{v^{2}} = \frac{\lbrack {{f/2} + {ɛ\; \sin^{2}\overset{\_}{\theta}}} \rbrack^{2} - \{ {\frac{v^{2}}{v_{p_{0}}^{2}} - \lbrack {1 + {ɛ\; \sin^{2}\overset{\_}{\theta}}} \rbrack + {f/2}} \}^{2}}{2f\; {v^{2}( {ɛ - \delta} )}\sin^{2}\overset{\_}{\theta}}} \\{f = {1 - ( \frac{v_{s_{0}}}{v_{p_{0}}} )^{2}}}\end{matrix} $

where phase velocity v has roots of two magnitudes: one for quasi-Pwaves, and the other for quasi-SV waves θ is the angle between thewavefront normal and the axis of symmetry, and other parameters aredefined in eq. [1].

According to the geometry shown in FIG. 13, the wavefront normal ({rightarrow over (n)}) and the axis of symmetry ({right arrow over (t)}) andthe angle in between take the following form:

$\quad\{ \begin{matrix}{\overset{arrow}{n} = {{\sin \; \theta \; \cos \; \varphi \overset{arrow}{i}} + {\sin \; \theta \; \sin \; \varphi \overset{arrow}{j}} + {\cos \; \theta \overset{arrow}{k}}}} \\{\overset{arrow}{t} = {{\sin \; \theta_{0}\cos \; \varphi_{0}\overset{arrow}{i}} + {\sin \; \theta_{0}\sin \; \varphi_{0}\overset{arrow}{j}} + {\cos \; \theta_{0}\overset{arrow}{k}}}} \\\begin{matrix}{{\cos \overset{\_}{\theta}} = \frac{\overset{arrow}{n} \cdot \overset{arrow}{t}}{{n}{t}}} \\{= {{\sin \; \theta \; \cos \; \varphi \; \sin \; \theta_{0}\cos \; \varphi_{0}} + {\sin \; \theta \; \sin \; \varphi \; \sin \; \theta_{0}\sin \; \varphi_{0}} + {\cos \; \theta \; \cos \; \theta_{0}}}}\end{matrix}\end{matrix} $

where θ₀ is the tilt of the axis of symmetry with respect to thevertical in a TI medium, and φ₀ is the azimuth of the axis of symmetry.Recognizing that:

$\quad\{ \begin{matrix}{{\sin \; \theta \; \cos \; \varphi} = {k_{x}{v/\omega}}} \\{{\sin \; \theta \; \sin \; \varphi} = {k_{y}{v/\omega}}} \\{{{\cos \; \theta} = {k_{z}{v/\omega}}},}\end{matrix} $

the following fourth-order dispersion relations can be derived:

ω⁴ −v _(p) ₀ ²[(1+2ε)f ₂ +f ₁−(f−1)f ₃]ω² +v _(p) ₀ ⁴[2f(ε−δ)f ₁ f₂−(f−1)((1+2ε)f ₂ +f ₁)f ₃]=0  [27]

6. Establishing a Fourth-Order PDE in TTI Media

Multiplying both sides of the above fourth-order dispersion relationwith a scalar wavefield P, and converting the frequency-wavenumberoperators to the time-space domain, the fourth-order PDE for TTI/VTImedia takes of the form of eq. [8].

7. Establishing a Second-Order 2×2 System of PDEs for TTI Media

The above fourth-order pseudo-acoustic PDE for TTI media can be solvedby the 2×2 time- and space-domain PDE system by eq. [3]:

${Where}\mspace{14mu} \{ \begin{matrix}{v_{h} = {v_{p_{0}}\sqrt{1 + {2ɛ}}}} \\{{v_{n} = {v_{p_{0}}\sqrt{1 + {2\delta}}}},}\end{matrix} $

the 2×2 system of PDEs can also take an equivalent form in terms ofhorizontal velocity v_(h) and normal-moveout (NMO) velocity v_(n).

In a 2D medium as a special case, the above PDEs still remain valid withthe following simplified spatial derivative operators:

$\quad\{ \begin{matrix}{f_{1} = ( {{\sin^{2}\theta_{0}\frac{\partial^{2}}{\partial x^{2}}} + {\cos^{2}\theta_{0}\frac{\partial^{2}}{\partial z^{2}}} + {\sin \; 2\theta_{0}\frac{\partial}{\partial x}\frac{\partial}{\partial z}}} )} \\{f_{2} = ( {{\cos^{2}\theta_{0}\frac{\partial^{2}}{\partial x^{2}}} + {\sin^{2}\theta_{0}\frac{\partial^{2}}{\partial z^{2}}} + {\sin \; 2\theta_{0}\frac{\partial}{\partial x}\frac{\partial}{\partial z}}} )}\end{matrix} $

8. Establishing a Second-Order 3×3 System of PDEs for TTI Media

As an alternative, the fourth-order pseudo-acoustic PDE for TTI mediacan also be solved by a 3×3 time- and space-domain PDEs in eq. [4] orits equivalents using a different linear combination.

Embodiments of the present invention can be implemented on eitherco-processor accelerated architectures, such asField-Programmable-Gate-Arrays (FPGAs), Graphics-Processing-Units(GPUs), Cells, or general-purpose computers. The present inventionprovides apparatus and general-purpose computers and/or co-processorsprogrammed with instructions to perform a method for the presentinvention, as well as computer-readable media encoding instructions toperform a method of the present invention. A system for performing anembodiment of the present invention is schematically illustrated in FIG.14. A system 52 includes a data storage device or memory 54. The storeddata may be made available to a processor 56, such as a programmablegeneral purpose computer. The processor 56 may include interfacecomponents such as a display 58 and a graphical user interface 60. Thegraphical user interface (GUI) may be used both to display data andprocessed data products and to allow the user to select among optionsfor implementing aspects of the method. Data may be transferred to thesystem 52 via a bus 62 either directly from a data acquisition device,or from an intermediate storage or processing facility (not shown).

Although the invention has been described in detail for the purpose ofillustration based on what is currently considered to be the mostpractical and preferred embodiments, it is to be understood that suchdetail is solely for that purpose and that the invention is not limitedto the disclosed embodiments, but, on the contrary, is intended to covermodifications and equivalent arrangements that are within the spirit andscope of the appended claims. For example, though reference is madeherein to a computer, this may include a general purpose computer, apurpose-built computer, an ASIC programmed to execute the methods, acomputer array or network, or other appropriate computing device. As afurther example, it is to be understood that the present inventioncontemplates that, to the extent possible, one or more features of anyembodiment can be combined with one or more features of any otherembodiment.

1. A computer-implemented method for pseudo acoustic quasi-P wavepropagation which remain stable in variable tilt anisotropic media andis not limited to weak anisotropic conditions, the method includes;acquiring a seismic exploration volume for a subsurface region ofinterest; determining a modeling geometry for the seismic explorationvolume; propagating at least one wavefield through the seismicexploration volume utilizing the modeling geometry for initialconditions and preventing the accumulation of energy along the axis ofsymmetry of anisotropic regions within the seismic exploration volumeand ensuring positive stiffness coefficients in the stress-strainrelations utilizing finite quasi-S wave velocities thereby producing astable wavefield; and utilizing the stable wavefield to generatesubsurface images of the subsurface region of interest.
 2. The method ofclaim 1 wherein propagating at least one wavefield through the seismicexploration volume includes an artificial quasi-shear wave velocitybeing greater or equal to zero along the axis of symmetry of anisotropicregions within the seismic exploration volume.
 3. The method of claim 1wherein propagating at least one wavefield through the seismicexploration volume includes restricting accumulation of quasi-shearwaves along the axis of symmetry for the seismic exploration volume. 4.The method of claim 1 wherein a plurality of wavefields are propagatedthrough the seismic exploration volume.
 5. The method of claim 1 whereinpropagating at least one wavefield through the seismic explorationvolume includes utilizing one of reverse time migration, a wave-equationbased migration, Gaussian beam migration or Kirchhoff migration.
 6. Themethod of claim 1 which includes propagating wavefields forwards andbackwards through the seismic exploration volume and applying imagingconditions to the forward and backward wavefields or equivalent Green'sfunctions to derive the subsurface images.
 7. The method of claim 6wherein the step of applying imaging conditions includes crosscorrelation between the forward and backward wavefields or equivalentGreen's functions to derive the subsurface images.
 8. The method ofclaim 7 wherein the step of applying imaging conditions includes atleast one of illumination normalization, reflection-angle domain gathergeneration and phase-amplitude compensation.
 9. The method of claim 1wherein the imaging output geometry includes common offset, commonazimuth and common reflection-angle domains.
 10. The method of claim 1wherein propagating at least one wavefield includes utilizing at leastone of delayed shot, plane-wave and phase encoding.
 11. The method ofclaim 1 wherein propagating at least one wavefield includes utilizing atleast one of normal moveout velocity, horizontal velocity and Thomsenparameters.
 12. A system configured to perform pseudo acoustic quasi-Pwave propagation which remain stable in variable tilt anisotropic mediaand is not limited to weak anisotropic conditions, the systemcomprising: a data storage device having computer readable dataincluding a seismic exploration volume for a subsurface region ofinterest, a processor, configured and arranged to execute machineexecutable instructions stored in a processor accessible memory forperforming a method comprising: determining a modeling geometry for theseismic exploration volume; propagating at least one wavefield throughthe seismic exploration volume utilizing the modeling geometry forinitial conditions and preventing the accumulation of energy along theaxis of symmetry of anisotropic regions within the seismic explorationvolume and ensuring positive stiffness coefficients in the stress-strainrelations utilizing finite quasi-S wave velocities thereby producing astable wavefield; and utilizing the stable wavefield to generatesubsurface images of the subsurface region of interest.
 13. The systemof claim 12 wherein propagating at least one wavefield through theseismic exploration volume includes an artificial quasi-shear wavevelocity being greater or equal to zero along the axis of symmetry ofanisotropic regions within the seismic exploration volume.
 14. Thesystem of claim 12 wherein propagating at least one wavefield throughthe seismic exploration volume includes restricting accumulation ofquasi-shear waves along the axis of symmetry for the seismic explorationvolume.
 15. The system of claim 1 wherein a plurality of wavefields arepropagated through the seismic exploration volume.